home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
newmat03.lha
/
newmat03
/
cholesky.cxx
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-08
|
916b
|
38 lines
//$$ cholesky.cxx cholesky decomposition
// Copyright (C) 1991: R B Davies and DSIR
#define WANT_MATH
#include "include.hxx"
#include "newmat.hxx"
/********* Cholesky decomposition of a positive definite matrix *************/
// Suppose S is symmetrix and positive definite. Then there exists a unique
// lower triangular matrix L such that L L.t() = S;
ReturnMatrix Cholesky(const SymmetricMatrix& S)
{
int nr = S.Nrows();
LowerTriangularMatrix T(nr);
real* s = S.Store(); real* t = T.Store(); real* tk = t;
for (int i=0; i<nr; i++)
{
real* tj = t; real* ti = tk;
for (int j=0; j<=i; j++)
{
tk = ti; real sum = 0.0; int k = j;
while (k--) { sum += *tj++ * *tk++; }
if (i==j) *tk++ = sqrt(*s++ - sum);
else *tk++ = (*s++ - sum) / *tj++;
}
}
T.Release(); // not wanted if routine avoids return init
return T;
}